Import Libraries

library(raster)
## Loading required package: sp
library(rhdf5)
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-3

Import Canopy Height Model

# read LiDAR canopy height model
chm <- raster("../NEONdata/D17-California/TEAK/2013/lidar/TEAK_lidarCHM.tif")

View CHM

# assign chm values of 0 to NA--take all of the 0 pixels and make them NA instead--so you can analyze all trees, for example
chm[chm==0] <- NA

# do the values in the data look reasonable?
plot(chm,
     main="Canopy Height \n LowerTeakettle, California")

image(chm,
      main = "These are just pixels and will stretch the space")

hist(chm,
     xlab="Tree Height")

Import Aspect Data

aspect <- raster("../NEONdata/D17-California/TEAK/2013/lidar/TEAK_lidarAspect.tif")
plot(aspect,
     main="Aspect for Teakettle Field Site")

## Classify North and South Facing Slopes

Create Classification Mstrix

# first create a matrix of values that represent the classification ranges
# North facing = 1
# South facing = 2
class.m <- c(0, 45, 1,
             45, 135, NA,
             135, 225, 2,  
             225 , 315, NA,
             315, 360, 1)
class.m
##  [1]   0  45   1  45 135  NA 135 225   2 225 315  NA 315 360   1
# reshape the object into a matrix with columns and rows
rcl.m <- matrix(class.m, 
                ncol=3, 
                byrow=TRUE)
rcl.m
##      [,1] [,2] [,3]
## [1,]    0   45    1
## [2,]   45  135   NA
## [3,]  135  225    2
## [4,]  225  315   NA
## [5,]  315  360    1

Reclassify the Raster

# reclassify the raster using the reclass object - rcl.m
asp.ns <- reclassify(aspect, 
                     rcl.m)
# plot outside of the plot region
# make room for a legend
par(xpd = FALSE, mar=c(5.1, 4.1, 4.1, 4.5))
# plot
plot(asp.ns,
     col=c("white","blue","orange"), # hard code colors, unclassified (0)=white,
         #N (1) =blue, S(2)=green
     main="North and South Facing Slopes \nLower Teakettle",
     legend=F)
# allow legend to plot outside of bounds
par(xpd=TRUE)
# create the legend
legend((par()$usr[2] + 20), 4103300,  # set x,y legend location
       legend = c("North", "South"),  # make sure the order matches the colors, next
       fill = c("blue", "orange"),
       bty="n") # turn off border

## Export a Geotiff

# export geotiff
writeRaster(asp.ns,
            filename="../NEONdata/outputs/TEAK/Teak_nsAspect2.tif",
            format="GTiff",
            ##compress file
            options="COMPRESS=LZW",
            ##remote sensing standard for NA
            NAflag = -9999)

Create a Raster Mask–this creates a mask layer, where we have the ndvi data only for pixels that have a north or south facing slope

ndvi <- raster("../NEONdata/D17-California/TEAK/2013/spectrometer/veg_index/TEAK_NDVI.tif")
plot(ndvi,
     main = "NDVI for Teakettle field Site")

# mask data
nFacing.ndvi <- mask(ndvi, 
                     asp.ns)
plot(nFacing.ndvi)